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ABSTRACT 

Hydrodynamical calculations in three space dimensions of the collapse of an isothermal, 
rotating 1 Mq protostellar cloud are presented. The initial density stratihcation is a 
power law with density p oc with p = 1. The case of the singular isothermal sphere 
{p = 2) is not considered; however p = 1 has been shown observationally to be a good 
representation of the density distribution in molecular cloud cores just before the beginning 
of collapse. The collapse is studied with two independent numerical methods, an SPH code 
with 200,000 particles, and a hnite-difference code with nested grids which give high spatial 
resolution in the inner regions. Although previous numerical studies have indicated that 
such a power-law distribution would not result in fragmentation into a binary system, both 
codes show, in contrast, that multiple fragmentation does occur in the central regions of 
the protostar. Thus the process of binary formation by fragmentation is shown to be 
consistent with the fact that a large fraction of young stars are observed to be in binary 
or multiple systems. 
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1 INTRODUCTION 


This paper considers a long-standing problem regarding binary formation by fragmentation 
of a collapsing, rotating clond, namely, that observed molecnlar clond cores in star forming 
regions are snfhciently centrally condensed so that they wonld not fragment according to 
n um erical calcnlations, yet these regions have a high proportion of binary systems. Fnller 
& Myers (1992), in a snrvey of clond cores, some of which are in Tanrns, fonnd a density 
distribntion p oc r“^, where p ~ 1.6, independent of whether the cores had IR sonrces in 
them or not. The snrvey by Ward-Thompson et al (1994), in the dnst continnnm, of cores 
withont infrared sonrces showed p ~ 1.25 in the inner regions and p ~ 2.0 in the onter 
regions, where the ‘inner’ region refers to r < 0.05 pc. One of these cores, L1689B, has been 
reobserved at 1.3 mm at higher spatial and angnlar resolntion (Andre, Ward-Thompson, 
& Motte 1996) and the resnlt has been conhrmed: the central regions have p ^ 1. In 
addition, theoretical stndies (Lizano & Shn 1989; Tomisaka, Ikenchi, & Nakamnra 1990; 
Ciolek & Monschovias 1994; Basn & Monschovias 1994) of the qnasistatic evolntion of 
magnetically snpported molecnlar clonds end np with cores having a centrally condensed 
distribntion at the onset of collapse. 

However, according to previons three-dimensional hydrodynamical calcnlations of the 
collapse of snch cores (Boss 1987; Myhill & Kanla 1992) they do not fragment if they are 
initially nniformly rotating and have p = 1. The same conclnsion is reached for clonds 
with p = 2, an often-nsed initial condition for n um erical calcnlations of star formation, and 
in this case the conclnsion is bolstered by a linear stability analysis (Tsai & Bertschinger 
1989). Althongh Myhill & Kanla (1992) do hnd fragmentation with an initially differen¬ 
tially rotating clond (p = 1 or 2), there is no observational evidence that clonds at the 
onset of collapse are differentially rotating, and fnrthermore, it is expected from theory 
(Basn & Monschovias 1994) that the initial clonds will be in near-nniform rotation becanse 
of the effects of magnetic braking. 

From what has jnst been said one wonld dednce, for example for the Tanrns-Anrigae 
star-forming region, that very few of the formed stars wonld be in binary systems. However 
in fact almost all of the stars in that region are mnltiple, as dednced from the observations 
of several gronps (Ghez, Nengebaner, & Matthews 1993; Leinert et al 1993; Simon et 
al 1995; Brandner et al 1996; Mathien 1994). In other star-forming regions (Brandner 
et al 1996) the binary fraction is not qnite so high as in Tanrns, bnt nevertheless it is 
comparable to that among nearby solar-type main-seqnence stars (Dnqnennoy & Mayor 
1991). It is a generally accepted conclnsion that most of these binaries mnst have formed 
dnring the protostellar collapse, perhaps by fragmentation; Fnller, Ladd, & Hodapp (1996) 
have recently discovered a binary protostar in the earliest stages of formation. 

In the present paper we consider the qnestion of whether fragmentation dnring pro¬ 
tostellar collapse is, in fact, a viable mechanism for binary formation. Althongh many 
previons stndies of three-dimensional protostar collapse have resnlted in fragmentation 
(Boss 1993, 1996; Bnrkert & Bodenheimer 1993,1996, Bonnell & Bastien 1993, Bonnell & 
Bate 1994, Monaghan 1994; Miyama 1992; Tnrner et al 1995; Nelson & Papaloizon 1993; 
Sigalotti & Klapp 1994), the initial assnmed density variation with distance from the centre 
was generally either nniform or a mild exponential (Boss 1991, 1993),that is, less centrally 
condensed than the distribntions snggested by observations. We calcnlate the collapse of 
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isothermal, rotating clouds, without magnetic helds, with two different three-dimensional 
hydrodynamic codes, starting with the density distribution p = 1 and with uniform rota¬ 
tion. The initial density distribution is assumed to be spherically symmetric with a small 
nonaxisymmetric perturbation. The purpose of our calculation is to calculate the collapse 
with much higher spatial resolution and to carry it over a longer time period than earlier 
numerical calculations to see if fragmentation does in fact occur. We also test to see if 
fragmentation properties depend on the assumed angular moment um of the cloud. 


2 THE COMPUTATIONAL METHOD AND INITIAL CONDITIONS 
2.1 Grid code 

The finite difference method is essentially the same as that described by Burkert & Bo- 
denheimer (1993, 1996). The calculations are performed on a 3-dimensional Eulerian, 
Cartesian grid; the advection scheme is based on second-order monotonic transport (van 
Leer 1977). The full computational region is represented by a standard grid, composed of 
64^ grid cells equally spaced in all directions. For improved resolution of the inner regions, 
four Cartesian nested concentric subgrids were superimposed on the standard grid, giving 
a ratio of total cloud size to the size of the smallest zone of 1000 or 2000. The linear scale 
on a given subgrid is reduced by a factor of 2 or 4 with respect to the next larger grid. 
The grid structure is set up at the beginning of the calculation and left fixed during the 
entire run. An artificial viscosity of the type described by Colella & Woodward (1984) is 
added; its main effect is to suppress a rapid increase in the density at centres of fragments. 
Unlike the calculations of Burkert & Bodenheimer (1996) there is no symmetry assumed 
with respect to the z-axis, although the configuration is assumed to be symmetric with 
respect to the equatorial plane. 

An important difference between the calculations reported in this paper and most 
previous fragmentation calculations is that an additional numerical stability criterion is 
incorporated. Klein, Truelove, & McKee (1996) have pointed out that numerical fragmen¬ 
tation could occur in a calculation if the local zone size is larger than the local Jeans length. 
In their adaptive mesh calculation they simply dynamically rezone to a sufficiently high 
level of accuracy so that the criterion is always satisfied. In our calculations, to avoid using 
an excessively large number of zones, we introduce an artificial heating at high densities, 
so that if the mass of a zone exceeds the local Jeans mass, the zone is heated until its mass 
is less than the local Jeans mass. The heating procedure affects only the regions of highest 
density, which in general are already part of fragments, and tests of calculations with and 
without heating show very little difference in the results. 

The Jeans mass is given by (Spitzer 1978) as p/j where Ij is the local Jeans length; 

thus 

Mj = 2.44 X g (1) 

p 

If the zone size is A, then the requirement that the zone mass pA^ < Mj defines a critical 
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density pc- 


( 2 ) 


3.905 X _o 

P<Pc^ -“ S cm 

If p > Pc, then the zone is heated so that its new sonnd speed Cg is 

. 2.561 X 2 -2 

ct = - - cm s 

Vp 


(3) 


Here Rg is the gas constant and rjp is a safety factor, taken to be 0.5, so that the heating 
starts at somewhat lower densities than pc. 

2.2 SPH code 

The SPH calcnlations were performed nsing a three-dimensional code based on a version 
originally developed by Benz (1990; Benz et al. 1990). The standard form of artihcial 
viscosity is nsed (Monaghan & Gingold 1983; Benz 1990; Monaghan 1992), with the pa¬ 
rameters Ov = 1 and /?v = 2. The smoothing lengths of particles are variable in time 
and space, snbject to the constraint that the nnmber of neighbonrs for each particle mnst 
remain approximately constant at iVneigh = 50. The SPH eqnations are integrated nsing a 
second-order Rnnge-Kntta-Fehlberg integrator. 

Two major modihcations have been made to the original code. First, individnal time 
steps are nsed for each particle (Bate, Bonnell & Price 1995) rather than a single time step 
for all particles. The resnlt is a great saving in compntational time for simnlations where 
there is a large density contrast (e.g. fragmentation calcnlations). Second, gravitational 
forces and a particle’s nearest neighbonrs are fonnd by nsing either a tree, as in the original 
code, or the special-pnrpose GRAvity-piPE (GRAPE) hardware. The implementation of 
SPH nsing the GRAPE closely follows that described by Steinmetz (1996). Using the 
GRAPE attached to a Snn SPARGstation 20 workstation typically resnlts in a factor of 
5 improvement in speed over the workstation alone. The calcnlation presented here was 
performed nsing the GRAPE. 

We nse the polytropic eqnation of state P = Kp'^, where AT is a constant that depends 
on the entropy of the gas. The polytropic constant 7 varies with density as 


7 = 1, P < Pc 2 = 10 g cm 

5 

7 = 3 , P> Pc2 = 10“^^ g cm“^ (4) 

and K is dehned snch that, when the gas is isothermal {p < Pc 2 ), K = with the 
sonnd speed Cg = 2.037 x 10^ cm s“^. Heating the gas when p > pc 2 serves two pnrposes. 
First, with a pnrely isothermal eqnation of state, a non-rotating fragment collapses to 
inhnite density. Thns, a fragmentation calcnlation mnst be stopped when the hrst fragment 
is formed. If the gas is heated, however, the collapse of a fragment is halted and the 
calcnlation can be followed fnrther. Second, heating of the gas is reqnired so that the 
minimnm mass that can be resolved (~ 2 A^neigh times the particle mass) is always less 
than a Jeans mass. If this criterion is not met fragmentation may be artihcially indnced 
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(Bate & Burkert, in preparation). This criterion is similar to that nsed for the grid 
code. For the calcnlation presented in this paper 200,000 eqnal-mass particles are nsed. 
Thns, from eqnation 1, a Jeans mass is less than the mass of 2A^neigh particles when 
p>7 X 10“^^ g cm“^. Heating the gas for p > pc 2 ensnres a Jeans mass is always resolved. 

2.3 Initial conditions 

The protostar in the present calcnlation is assnmed to be isothermal in space and time. 
The initial conditions are specihed by the ratios a and /3, which are, respectively, the 
thermal and rotational energies divided by the absolnte valne of the gravitational energy, 
and by the angnlar momentnm distribntion, the size and form of the initial pertnrbation, 
and the form of the density distribntion. The latter is taken to be 

with p = 1. Here r is the distance to the origin and po and tq are constants. The form of 
the pertnrbed prohle is 

Pvi'^) = Pi^) [1 + cos(2(/))] (6) 

where ai is the amplitnde of the pertnrbation and cf) is the azimnthal angle abont the 
rotation (z) axis. In all cases presented here the ai is set to 0.1, the clond mass Mtot is 
set to 1.0 Mq, ropo = 0.1273 g cm“^, the angnlar velocity O is nniform, and the radins R 
of the original sphere is 5 xlO^® cm. The sonnd speed Cg = 2.037 x 10^ cm s“^; thns the 
valne of a for all calcnlations was 0.35. 


3 RESULTS 

3.1 Analytic discussion 

A spherically symmetric, isothermal gas clond with a power-law density distribntion as 
given by eqnation 5 is a somewhat nnphysical initial condition for 0 < p < 2 becanse there 
exists a critical radins Tcrit > 0 inside which the gas will initially expand even if centrifngal 
forces are neglected. Within Vcrit the pressnre gradient initially exceeds the gravitational 
force: 


^ GM{r) 
p dr ~ 

where M(r) is the total mass inside r and eqnality holds for r = Tcrit- Given the density 
distribntion of eqnation 5 and nsing eqnation 7, we hnd 


^ crit R 


/ 2p(3 — p)a 
V 3(5 - 2p) 


2-p 


( 8 ) 
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which is valid for 0 < p < 2. 

Figure 1 shows Vcrit/R and the fractional mass of the expanding region M{rcrit)/Mtot 
as a function of p for different values of a. Note that all power-law prohles with 0 < p < 2, 
even with small values of a, have an initially expanding core. For the present case, where 
a = 0.35 and p = 1, Vcrit/R = 0.156 and M(rcrit)/Mtot = 0.024. 


3.2 Fragmentation, grid code, high beta 

The hrst run was calculated with 13 = 0.23 on four subgrids, each 128 xl28 x 64 zones 
in the {x,y,z) directions, respectively, which had radii of 2.50 xlO^® cm, 1.25 xlO^® cm, 
6.25 xlO^® cm, and 3.125 xlO^^ cm. The resolution on the hnest grid is 4.88 xlO^^ cm. 
Although the initial central density was 6 xl0“^® g cm“^ the initial expansion phase, 
which affects the inner 25% of the radius, results in a reduction of the central density to 
about 2 xl0“^^ g cm“^. The density and velocity distributions on the hrst subgrid at a 
time 4.3 xlO^^ s are shown in Figure 2. At this time the expanding central region has 
nearly reached its maximum extent, and is about to resume collapse. The density prohle 
is markedly less steep in this region than in the original cloud. The free-fall time in the 
central regions is now 4.4 xlO^^ s. 

The collapse then proceeds through more than six orders of magnitude increase in 
the central density before any signs of fragmentation begin to appear. Figure 3 shows 
the development of the central region (the innermost subgrid) starting at 1.176 xlO^^ s 
from the beginning of the calculation. In Figure 3a a dehnite central density maximum is 
present. A small region around it has formed a battened disk structure and a spiral-arm 
pattern. At the end points of the spiral pattern one can already see the development of two 
condensations which accrete material until they become self-gravitating. Figure 3b shows 
the original central density maximum connected by inner spiral arms to the newly-formed 
binary. Outside this region an outer spiral pattern has developed and a circum-triple 
disk is beginning to form. This high-density disk has become well-developed a short time 
later (Fig. 3c). At that time the symmetry of the inner system is broken as a result of 
small numerical perturbations and the triple system rearranges itself into a more stable 
hierarchical system, with a close binary and a third object further out. At the outer edge 
of the disk two new density maxima are visible, triggered by the interaction of the outer 
spiral arms. One of these maxima hnally develops into an outer fragment, the other one 
is tidally disrupted and subsequently accreted by the outwards moving fragment of the 
previous inner triple system (Fig. 3d). At t = 1.217 x 10^^ s the masses of the inner binary 
components are 0.055 and 0.025 Mq. The mass of the third component of the initial triple 
system is 0.017 Mq. The fourth fragment has a mass of 0.011 Mq. 

At this time we stop the calculation; on the inner grid (scale 3 xlO^^ cm) 90% of the 
total mass is included in the four fragments. However only a small fraction (12.4%) of Mtot 
is included in this region. We expect further interactions of these fragments to take place 
and possibly further fragmentation to occur on larger scales. But at this point it is clearly 
established that fragmentation has occurred, that some of the fragments will survive, and 
that the end result will be a multiple system. 
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3.3 Fragmentation, SPH code, high beta 

To test whether or not the fragmentation is dependent on the nnmerical method, we 
perform the same calcnlation as in Section 3.2, bnt nse Smoothed Particle Hydrodynamics 
(SPH). We nse 200,000 eqnal-mass particles in order to resolve a Jeans mass at densities 
np to 10“^^ g cm“^ while maintaining an isothermal eqnation of state. Beyond 10“^^ g 
cm“^ the gas is heated so that a Jeans mass is always resolved. 

The early evolntion is identical to that calcnlated by the grid code, with an expansion 
in the centre of the clond. The initial resolntion in the centre is not as good as with the 
grid code; the maximnm density is only 5 x 10“^® g cm“^. This resnlts from the fact that, 
with SPH, the resolntion depends on the density. The smoothing length h, which gives 
the resolntion, is dehned so that 2h contains A^neigh particles. In the centre of the clond, 
althongh the density formally goes to inhnity as the radins goes to zero, the mass fraction 
is very small, and hence h mnst be large in order to contain A^neigh particles. Other than 
the difference in initial central density, however, the expansion is identical to that given by 
the grid code, with the same density minimnm of 2 x 10“^^ g cm“^ being reached before 
the expansion stops and the collapse begins. 

Dnring the collapse, the central density increases by more than six orders of magnitnde, 
in agreement with the grid code. A central density maxim um forms, and a disk of gas begins 
to form aronnd it. The evolntion of the central density maximnm and its disk is shown in 
Fignre 4. As soon as the disk appears aronnd the central density maximnm it is threaded by 
spiral arms (Fig. 4a). The arms continnally wrap-np, interact and reform becanse of their 
differential rotation. Ultimately, density maxima form at the ends of the two arms leading 
to a binary aronnd the central object (Fig. 4b). As higher angnlar momentnm material 
falls in, spiral arms extend beyond the binary into a circnm-triple disk. The symmetry 
of the triple system rapidly breaks, cansing the onter spiral arms to become asymmetric. 
Again, density maxima form at the ends of these arms, however, dne to the asymmetry, 
only one of these forms a fragment immediately; the other is temporarily disrnpted (Fig. 
4c). Eventnally, this disrnpted hfth density maxima does gather enongh material to form 
a hfth fragment (Fig. 4d). Note that a hfth fragment appears also temporarily in the grid 
code. There, however, the somewhat diherent orbital evolntion of the inner triple system 
leads to a merger (Fig. 3d). 

The collapse and fragmentation follows the same pattern as with the grid code: a 
central object is formed, snrronnded by a disk with spiral arms; fragments form at the end 
of the arms; the initial symmetry of the triple system breaks with one fragment moving 
ontwards and the other two components forming a close binary system; a circnm-triple 
disc forms, again threaded by spiral arms which fragment, with the details depending 
on symmetry breaking of the triple. The fact that this pattern of fragmentation has 
been produced with two completely different hydrodynamic codes demonstrates that the 
fragmentation is physical and not a numerical effect. Some details of the fragmentation 
differ between the two codes. The formation of the triple system occurs slightly earlier 
with the SPH code {t = 1.155 x 10^^ s) than the grid code {t = 1.176 x 10^^ s), before 
the disc has reached as large an extent. This earlier fragmentation is probably due to 
the different resolutions of the two codes, and the higher numerical noise that is present 
in the SPH code. The resolution of the two codes are almost identical at densities of 



1.0 X 10“^^ g cm“^ {h = 4.0 x 10^^ cm, while the smallest grid cell has A = 4.9 x 10^^ cm). 
However, because the resolution of SPH depends on the density, for p > 10“^^ g cm“^ 
(i.e. in the spiral arms) the resolution is actually higher with SPH. The higher numerical 
noise of the SPH code results from the use of the GRAPE hardware which has less than 
single precision, whereas the grid code calculation (using a CRAY) performs gravitational 
calculations with double precision. This difference in numerical noise also leads to the 
symmetry of the triple system being broken quicker with the SPH code than the grid code. 

We stop the calculation at the time 1.190 x 10^^ s (Fig. 4d). At this point the disk 
and fragments contain 12.4% of Mtot- The fragment masses are 0.050 and 0.020 Mq for 
the close binary and 0.018 Mq for the third component of the initial triple system. These 
values are in excellent agreement with the results of the grid code. The two fragments 
which formed later on the far-left and far-right of Figure 4d have masses of 0.017 Mq and 
0.003 Mq, respectively. As with the grid code, mergers between fragments and further 
fragmentation are probable as the remaining material falls in, but the survival of a multiple 
system is likely. 

3.4 Fragmentation, grid code, low beta 

The runs discussed in the previous sections were calculated with fairly large rotation {(3 — 
0.23). In order to check whether such a high value of [3 is crucial for fragmentation we 
have recalculated the same case with (3 reduced by a factor of 2. As the gas is now less 
rotationally supported we expect the system to become more centrally condensed. We 
therefore increased the resolution on the innermost grid by a factor of two, compared with 
the grid setup used in section 3.2. The four subgrids, each 128 xl28 x 64 zones in the 
{x,y,z) directions, respectively, now have radii of 2.50 xlO^® cm, 6.25 xlO^^ cm, 3.125 
xlO^^ cm, and 1.563 xlO^® cm. The resolution on the hnest grid is 2.44 xlO^^ cm. 

Figure 5 shows snapshots of the evolution, which is qualitatively the same as for the 
high-/? case, except that the scale of the region of initial fragmentation is reduced by a factor 
of 2. The initial development of spiral arms with density maxima at the ends is shown in 
Figure 5a. Figure 5b shows an inner triple system which formed in an analogous manner to 
the systems shown in Figures 3b and 4b. The symmetry then breaks as the original central 
density maximum approaches one of the new fragments, and a fourth fragment appears 
near the outer edge of the disk (Fig. 5c). A short time later hve fragments are present (Fig. 
5d) and the surrounding circum-fragment disk has expanded considerably. The calculation 
is stopped before the disk has grown to the size and mass where induced fragmentation 
on a larger scale (such as that shown in Fig. 3d) would be expected. However we obtain 
again multiple fragmentation; hve distinct fragments, including one of very low mass, have 
clearly been established (Fig. 5d.) The typical separations are now of order 10^® cm as 
compared with 4 x 10^^ cm at the end of the previous runs. The masses of the fragments 
in Figure 5d are (from left to right) .016, .017, .03, .002, and .016 Mq, slightly lower in 
total than in the previous case but including practically the entire mass on the innermost 
subgrid (scale 1.56 xlO^^ cm). Note also that in this case the timescale of fragmentation 
is slightly shorter. This case was not run further because at t = 1.007 x 10^^ s it is clear 
that fragmentation does occur even with the lower value of (3. 

The fragmentation is very likely not to be of numerical origin, since it all takes place 
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on the innermost grid and since heating is inclnded to keep zone masses below the Jeans 
mass. In fact, heating plays only a very minor role in the seqnence of events. From 
eqnation (2), the critical density above which heating takes place is pc = 3.27 x 10“^^ g 
cm“^. The densities at the centres of the fragments in Fignre 5d are typically 1.6 xl0“^^ 
g cm“^. and p > pc in only a few zones near the centres of the established fragments. 
The temperatnre increases to 98 K at the highest densities, which of conrse wonld tend to 
snppress fragmentation, bnt only in material which has already fragmented. The heating 
effect thns jnst slows the increase in the central densities. However the density of typical 
material in the vicinity of the fragments, which is the same as that from which the fragments 
originally condense, is abont 10“^^ g cm“^. There the heating has no effect. 


4 CONCLUSIONS 

The main pnrpose of this paper has been to demonstrate that mnltiple fragmentation 
occnrs even from an initial power-law density distribntion which is similar to what is 
observed in molecnlar clond cores. Onr resnlts differ from those of previons n um erical 
simnlations which, for initially nniform rotation, did not show evidence for fragmentation. 
However the initial expansion phase, which is an important factor as it leads to central 
expansion and to a core with a flat density distribntion, had previonsly not been resolved. 
Thns in fact the actnal initial conditions for collapse are nnlikely to have p oc r~P where 
0 < p < 2 bnt rather a distribntion in which the central regions have a flatter prohle. If 
snch a conhgnration has solid-body rotation its angnlar momentnm distribntion will be 
somewhat different from that nsed here, bnt it is jnst as likely to fragment, as has been 
demonstrated by the calcnlations, for example, of Boss (1993). 

In any case, high enongh spatial resolntion is reqnired to bring ont the detailed strnc- 
tnre of the central regions. In addition, the calcnlations have to be carried ont for a 
snfhcient time beyond that when a condensed central object has formed so that a snffi- 
cient amonnt of mass with higher angnlar momentnm has accnmnlated in the central disk 
region. Onr resnlts show that at least three to fonr fragments with nneqnal masses form, 
with the details depending on the valne of /? and the exact time of symmetry breaking 
which depends on the n um erical noise. 

It is promising that two entirely independent n um erical techniqnes lead basically to 
the same resnlt: the initial formation of a central triple system, snbseqnent symmetry 
breaking, and the formation of an onter disk with additional fragmentation on a larger 
scale. As the farther evolntion of the system is expected to be highly chaotic it will 
depend on the n um erical techniqne and the local resolntion which are not exactly identical 
in the two codes. It still has to be investigated by longer-term calcnlations to what extent 
the hnal stage will depend on all these details. Another important qnestion is whether 
mnltiple fragmentation will occnr in the extreme case of an p = 2 power law which does 
not lead to an initially expanding core. 
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FIGURE CAPTIONS 

Figure 1. The radius iuside of which a spherical power-law cloud iuitially expauds, 
is showu (left pauel), iu uuits of the total cloud radius, as a f un ction of the power-law 
index p, for various values of a, the initial ratio of thermal energy to the absolute value 
of the gravitational energy. The right panel shows the mass fraction inside iCcrit, as a 
function of the same two parameters. From top to bottom, the curves correspond to 
a = 0.5, 0.4, 0.3, 0.2, and 0.1. 

Figure 2. Contours of equal density in the equatorial (x, y) plane at z = 0 at a time of 4.3 
xlO^^ s in the high-beta grid-code calculation. The hrst subgrid is shown. The maximum 
density is log Pmax — —16.67 and the contour interval is A log p = 0.03. Velocity vectors 
are shown with length proportional to speed; the maximum velocity Vmax is 3.58 xlO^ cm 
s“^. The linear scale is given in cm. 

Figure 3. Evolution of the high-beta case with the grid code. The innermost grid is 
shown. Symbols and curves have the same meaning as in Fig. 2. (a; upper left) t = 1.176 
xlO^^ s; log Pmax — —10.3; A log p = 0.166; Vmax = 1-95 X 10^ cm s“^. (b; upper right) t 
= 1.190 xlO^^ s; log Pmax = —10.2; A log p = 0.17; Vmax = 2.57 x 10^ cm s“^. (c; lower 
left) t = 1.209 xlO^^ s; log pmax = —10.1; A log p = 0.173; Vmax = 3.17 x 10^ cm s“^. 
(d; lower right) t = 1.217 xlO^^ s; log pmax ~ —10.1; A log p = 0.173; Vmax — 3.67 x 10^ 
cm s“^. 

Figure 4. Evolution of the high-beta case with the SPH code. Symbols and curves have 
the same meaning as in Fig. 2. (a; upper left) t=1.135 x 10^2 s; log pmax = -10.0; A log 
p = 0.25; Vmax = 2.44 x 10^ cm s“^. (b; upper right) t=1.160 x 10^^ s; log pmax = —9.7; A 
log p = 0.25; Vmax = 2.88 x 10^ cm s“^. (c; lower left) t=1.175 x 10^^ s; log pmax = —9.6 A 
log p = 0.25; Vmax = 3.66 x 10^ cm s“^. (d; lower right) t=1.190 x 10^^ s; log pmax — —9.5; 
A log p = 0.25; Vmax — 3.38 x 10^ cm s“^. 

Figure 5. Evolution of the low-beta case with the grid code. The innermost grid is shown. 
Symbols and curves have the same meaning as in Fig. 2. (a; upper left) t = 9.637 xlO^^ s; 
log Pmax — —9.65; A log p = 0.219; Vmax — 1-98 X 10^ cm s“^. (b; upper right) t = 9.706 
xlO^^ s; log Pmax = —9.498; A log p = 0.224; Vmax = 1-88 x 10^ cm s“^. (c; lower left) t 
= 9.840 xlO^^ s; log Pmax = —9.46; A log p = 0.225; Vmax = 3.01 x 10^ cm s“^. (d; lower 
right) t = 1.007 xlO^^ s; log Pmax — —9.431; A log p = 0.226; Vmax = 3.41 x 10^ cm s“^. 
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